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Abstract Keywords Granular Medium • Drag 
Force • Event Driven Simulation with Friction 

We investigate the dynamics of an intruder pulled 
by a constant force in a dense two-dimensional granu- 
lar fluid by means of event-driven molecular dynamics 
simulations. In a first step, we show how a propagating 
momentum front develops and compactifies the system 
when reflected by the boundaries. To be closer to re- 
cent experiments [l][2]) we then add a frictional force 
acting on each particle, proportional to the particle's 
velocity. We show how to implement frictional motion 
in an event-driven simulation. This allows us to carry 
out extensive numerical simulations aiming at the de- 
pendence of the intruder's velocity on packing fraction 
and pulling force. We identify a linear relation for small 
and a nonlinear regime for high pulling forces and in- 
vestigate the dependence of these regimes on granular 
temperature. 
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1 Introduction 

The response of an intruder to a pulling force is a 
rather versatile tool to study the local nonequilibrium 
dynamics in various complex fluids, such as glasses [3], 
colloids [HIS] and granular media JLi2,iiX- The mea- 
sured force-velocity relations reveal striking nonlinear 
behaviour close to the glass and/or jamming transition. 
In the early experiments on colloids evidence was given 
for a threshold force even in the fluid regime, whereas 
mode-coupling theory predicts such behaviour only in 
the glassy phase. In the experiment of Ref. [TJ[5] two 
transitions are observed: the first one, called fluidiza- 
tion, separates a regime of continuous to intermittent 
motion of the intruder. It occurs below the jamming 
point (the second transition) and depends on the ap- 
plied pulling force. Experimentally it is observed even 
for the smallest pulling force with the possibility of a 
dynamic transition inherent to the vibrated granular 
fluid for zero applied force. In Ref. [Bl , the dynamics of 
an intruder was simulated near the jamming point. One 
result of this study are velocity-force relations which 
are linear for packing fraction, rj < 0.833 and become 
nonlinear for rj still closer to the jamming point. 

In contrast to [8] , we discuss a stochastically driven 
system, describing a fluidized granular state - similar to 
the experimental setup in [Tl[2] . A recent mode-coupling 
theory 131 predicts that such a "thermalized" granular 
fluid undergoes a glass transition at a packing fraction 
below the jamming point. This transition is different 
from both, the jamming transition at zero temperature 
and the glass transition for either Newtonian or Brow- 
nian dynamics. 

In this paper we use event driven simulations to an- 
alyze the dynamics of an intruder in a two-dimensional 
system of hard disks close to the glass transition. We 
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compute force-velocity correlations in the linear and 
nonlinear regime, extract the mobility for the linear 
regime and discuss scaling for the nonlinear regime. 
Moreover we discuss the dependence on the granular 
temperature. 

2 Model 

We consider a bidisperse system of A'' hard disks in an 
area A. Particles' positions and velocities are denoted 
by {r^lili and {vj.^^. The size ratio Rs/Rb = 4:5 
of small to big particles is chosen as in the experiments 
of Refs. [Hll]. The mass ratio follows from the disk like 
shape of the particles as ms/rrib = 16/25. The intruder 
with position Tq and velocity Vq is chosen twice as large 
as the small particles, Rq = 2Rs and mg = Amg. 

The particles collide inelastically so that in each 
collision energy is dissipated while momentum is con- 
served. The simplest model ignores rotational degrees 
of freedom, allowing for normal restitution only. The 
collision rules for two colliding particles (say i and 2) 
is given in terms of their relative velocity g := vi — V2 

(g-n)' = -e(g.n). (1) 

where e is the coefficient of restitution and n denotes 
the unit vector n := (ri — r2)/|ri — at contact. 

In between collisions the particles perform Brown- 
ian motion, due to friction with a surrounding medium 
or with the bottom plate. We model this by a frictional 
force in the equation of motion: Ffr — — 7mv. For all 
finite values of the friction coefficient 7, momentum is 
lost, whereas for 7 = 0, the total momentum is con- 
served. 

Inelastic collisions give rise to a loss of energy oc 
^—^T, where T denotes the granular temperature T = 

^f^i iTiivf . To balance this dissipation the entire 
system is driven stochastically (like an air-fluidized bed 
[TUlimiT^ or a vibrating bottom plate |131I141 [T5]) by 
instantaneous kicks. The momentum of the i-th particle 
Pi{t) is changed according to 

p'At) = P^{t) + PdrUt) (2) 

with (<ef (i)^f (^')) = 3,jSc,0S{t-t') and zero mean {^f{t)) 
0, and a = X, z denoting the cartesian components. 
In the frictional case, 7 ^ 0, we kick single particles, 
whereas for the frictionless case, 7 = 0, we kick neigh- 
bouring particles with equal and opposite momenta in 
order to conserve the total momentum locally [16]. 

We want to investigate the dynamics of a tagged 
particle under the action of a deterministic pulling force 
F. This so called intruder with position Tq and velocity 
Vq is subject to systematic kicks 

p'o{t)^Po{t)+FAt (3) 



mimicking a constant force acting on the intruder. The 
frequency of these systematic kicks {At)~^ is chosen 4 
orders of magnitude higher than the frequency of the 
stochastic driving and the collision frequency (see be- 
low). The systematic force does work on the system, 
injecting momentum and energy. 

Combining inelastic collisions, stochastic driving, pullin 
force and friction, we arrive at the following equation 
of motion 



dpl 
dt 



dt 



(4) 



coll 



We aim to prepare the granular fluid without pulling 
force to be under stationary conditions. This can be 
achieved by balancing energy dissipation due to fric- 
tion and inelastic collisions by the stochastic driving. 
Without external forcing the mean velocity of the par- 
ticles vanishes, so that the granular temperature, T = 
2^ X^iLi '^i'^'ii is just the average kinetic energy of the 
particles. Its time rate of change can be deduced from 
Eq. |4] by multiplying with pf (t) yielding 



1 d_ 

'2 dt 



ipf? 



1 d 
'2 di 



{Pff 



The explicit calculation of the time derivations on the 
rhs are beyond the scope of this paper but may be found 
elsewhere [TTlfTS] . The time rate of change of the gran- 
ular temperature then reads as follows: 



dr 1 

- = -27r-a;^- 



-T + /dr 



pI, 



2 2mcff 
where the effective mass is given by 



(6) 



simplicity, we choose the driving frequency /dr equal 
to the Enskog frequency oje [H] and measure length 
and mass in units, such that Rg = 1 and ~ 1. 
In the stationary state, ^ = 0, the amplitude of the 
stochastic driving is given by: 

-T. (7) 



Plr 



2jT 



1 



2mcff uje 2 
We prepare our system in a stationary state with T = 1 , 
yielding a numerical value for pdr depending on 7, oje 
and e. 



3 Simulations 

3.1 Implementation as event-driven simulation 

In order to apply an event-driven algorithm to the dy- 
namics, as described by Eq. (|4|), we need to gener- 
alize the code to account for damping (7 ^ 0). As 
usual, events include collisions of particles, wall colli- 
sions, subbox-wall collisions and (discrete) driving events 
(kicks). Standard event driven algorithms calculate the 
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time of an upcoming event and advance all particles to 
that time under the assumption that the particle mo- 
tion is ballistic. Hence, in between events, the velocities 
do not change. The condition for a collision is Ri + Rj — 
\ri{t') — rj{t*)\, i.e., the difference in the particles' tra- 
jectories at the collision time t' must be equal to the 
sum of their radii. Plugging in the trajectories subject 
to ballistic motion ri_j{t') = rij{to) + Vij(to)(t* — to), 
one finds a quadratic equation, that must have a posi- 
tive solution for (<* — to) if and only if a collision will 

occur uniin]- 

In our case, the velocities decrease as Vij(t*) = 
Vi,j(to) exp(— 7(t* — to)). Integration results in 

1 _ g-7(t*-to) 

r,,j{t*) = ri,j(tn)+Vij{to) . (8) 

7 

Inserting this into the condition for a collision, we find 
a similar condition for a collision except that (t* — to) 

is replaced by which is monotonically in- 
creasing in (t* — to). Since our (ballistic) event-driven 
code calculates (t* — to) we only have to replace the 
collision time by 

t*-to = --log(l-7.(f-to)). (9) 

7 

Hence, collisions will occur at the same place and in 
the same order but at a different time and with dif- 
ferent particles' velocities, provided the place of the 
collision is within reach of the particle. In detail, the 
maximum distance a particle is able to travel is limited 
because of the friction slowing the particle down. The 
maximal distance is given by r^°^ = hm(-._+oo ri(t*) = 
ri(to) + Vi(to)/7 > i"i(t*) that must be larger than the 
distance to the place at which the event will take place. 
Hence, (ri(t*) - r,(to)) < Vi(to)/7 or (f - to) < I/7. 
This is consistent with Eq. © which requires 
log (1 — 7 ■ (t* — to)) to be negative for a positive col- 
lision time. Advancing a particle to a different type of 
event ( e.g. kicks) is done in the same way. 

Changing an existing event-driven simulation in the 
way discussed enables us to simulate systems of hard 
disks and spheres subject to friction almost as fast as 
without friction. The systematic force on the intruder 
as described in Eq. [3] is implemented as frequent kicks 
on the intruder. In order to avoid an inelastic collapse, 
i.e., a infinite number of collisions in a finite time in- 
terval, we use the same method as described in [161 to 
circumvent it. 

3.2 Equilibration and data aquisition 

For packing fractions up to 77 = 0.8 we used the same 
data sets from [52] as initial conditions. For even denser 



systems, we used a compactified system acquired as de- 
scribed in section HH] 

The data shown in section 14.21 were obtained by 
first equilibrating 100 different configurations and then 
switching on the force on the intruder. The final ve- 
locity of the intruder was measured after stationarity 
was attained. The time to reach a stationary state de- 
pends on the packing fraction, the restitution coefficient 
and the applied force on the intruder. A typical trajec- 
tory and its corresponding fluctuating velocity for the 
largest packing fraction 77 = 0.8 and force F = 100000 
is shown the inset of Fig.|4l The intruder moves approx- 
imately 3i?o in a;-direction before fluctuating around its 
stationary average velocity. 

4 Results 

The time rate of change of the total momentum P = 
TO^Vi follows from Eq. (g]): 

P"(t)--7P"(t) + ^<5„,, (10) 
which is solved by 

P^{t) - P"(0)e-^* + Sa,.^{l - e-^*) (11) 

7iv 

Here we have used that the collisions as well as the 
random driving conserve momentum. In the frictional 
case with 7^0, the total momentum goes to a con- 
stant, Px ~ F/{Nj), whereas for the frictionless case 
with 7 = 0, the momentum grows linearly with time, 
Px — Ft / N . In the following we shall mainly discuss the 
first case, because the frictional model is closer to exper- 
iment. However, the frictionless case allows us to study 
the propagation of momentum in an inelastic fluid [231 

[21]. 

4.1 Frictionless state (7 = 0) 

To that end, we consider a system with aspect ratio 
^ = 5 : 1, = 5000 and fixed walls, which reflect the 
particles elastically. Momentum is conserved on aver- 
age except for the pulling force which constantly feeds 
(small) momenta into the system, which are propagated 
by collisions away from the source into all of the avail- 
able area. Below the jamming transition momentum 
transport is given by ballistic motion of particles as well 
as by collisions. If the intruder starts at the left hand 
side of the system and is pulled by the external force, 
then only particles in the neighbourhood of the intruder 
will feel the local momentum input for short times. The 
momentum given to the intruder is distributed among 
the particles in front (i.e. in the direction of the pulling 
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Fig. 1 Momentum transport with 7 = and hard wahs. 
A vector denoting the particle's velocity is assigned to each 
particle's position. (T = 0.1, £ = 0.9, F = 500, r) = 0.8, 
t = 2.6, 6.3, 8.5) 



i = 1.3125 
t = 3.8125 
t = 6.3125 




Fig. 2 Propagation of momentum wave; total momentum 
for 3 different times, in dimensionless units; (parameters as 
in Fig.[lll. 



force) of the intruder. Again, these particles distribute 
their momenta by colhsions as well as by ballistic trans- 
port. A front of particles carrying the momentum fed 
into the system by the intruder is formed (see Fig. [T] 
top), propagates through the system (see Fig. [1] center) 
and ultimately collides with the hard wall (see figure [T] 
bottom). At this instant the momentum is reflected by 
the wall and many particles end up in a highly com- 
pactified state with a packing fraction r] w 0.839. 

To analyze the momentum wave quantitatively, we 
plot Px averaged over z in Fig. [2] The propagation front 
is well defined, its center of mass moves with veloc- 
ity Vcm = 16.51 and it broadens with time such that 
it's width increases linearly with time, A — Vhit (with 
Vbr = 10.89). Both velocities increase with pulling force, 
while their ratio is approximately constant. Hence the 
observed propagation front cannot be identified with 
linear sound but is presumably a Shockwave in agree- 
ment with the observations of Ref . [23] . 



4.2 Force-velocity relation ( 7 7^ ) 



4-2.1 Linear Regime: Mobility of the intruder 

For a small pulling force we expect a linear relation 
between the velocity of the intruder and the force: 

VI = tiF, (12) 

defining the mobility fj, of the intruder. In Fig. ([3]) we 
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Fig. 3 Velocity versus pulling force for packing fractions 
0.3 < < 0.8, parameters chosen: 7 = l,e = 0.7 and as- 
pect ratio A = 1 : 2. 



In the following, we consider a system with N — 20000 
grains, subject to friction (7 = 1) with hard walls in the 
z-direction and periodic boundary conditions in the x- 
direction, allowing for a stationary current. We analyze 
the steady state motion of the intruder for small and 
large forces as a function of packing fraction for T = 1 
and subsequently compare these results to the same 
system with T = 0.04. 



plot the velocity of the intruder versus force and indeed 
do observe a linear regime for small pulling force. From 
the slope we extract the mobility which is shown in Fig. 

The breakdown of the mobility when the glass tran- 
sition is approached is clearly visible for both investi- 
gated e. The more inelastic system with e = 0.7 shows 
increased mobilities as compared to e = 0.9. Since the 
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Fig. 4 Mobility as a function of volume fraction for e = 0.9 
(stars) and e = 0.7 (squares). Inset: Typical time dependence 
of the X position and x velocity of the intruder for e = 0.9, 
F = 100000 and rj = 0.8. 

particles are more "sticky" , they tend to stay closer af- 
ter a collision than in the elastic limit, i.e., stronger den- 
sity fluctuations are inherent to systems with stronger 
inelasticity. This enlarges the accessible space for the 
intruder compared to the system with e = 0.9, increas- 
ing its average velocity. Moreover, the momentum of 
the intruder in the direction of the pulling force is re- 
duced due to collisions with other particles in front. 
In the center of mass frame, the intruder is reflected 
backwards. For the more inelastic system this kind of 
backscattering is less effective (see Eq. [T]), so that the 
velocity in the direction of the force and hence the mo- 
bility is larger for the more inelastic system. 

4-2.2 Nonlinear regime 

As the pulling force is increased, deviations from linear 
behaviour are expected and observed. To investigate 
these systematically we have applied forces in the range 
1 < F < 10^ and show our data in a scaling plot in Fig. 

Scaling velocities by v^i = Pdi/nT-cS and forces by 
^dr ■ V^^'^ collapses the data for large forces. The veloc- 
ity of the intruder scales algebraically with the pulling 
force, according to vj cx , (3 — 0.55. 

The crossover between the linear and the algebraic 
regime is accompanied by range of forces in which the 
intruder velocity increases superlinearly with the ap- 
plied force. This is observed at high packing fractions 
only, see Fig. [3] We conjecture that this shear thinning 
is due to the formation of vortices which can only be 
formed at high packing fractions, where momentum is 
conveyed almost instantaneously as the disks are very 
close to each other. For low packing fractions, the dis- 
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Fig. 5 Velocity-force relation for all pulling forces, e = 0.9. 

tance of neighbouring particles is too large for the for- 
mation of vortices because the damping 7 dissipates 
most of the momentum. 



4-2.3 Temperature dependence 

The crossover from linear to nonlinear response depends 
on the thermal velocity, uth — \/'^T/mo w 0.7. For 
small packing fractions, this crossover actually occurs 
at vth as can be seen e.g. in Fig. |2] for the smallest 
packing fractions. Decreasing the temperature is ex- 
pected to shrink the linear regime also for higher densi- 
ties. Hence we try to explore the force velocity relation 
in the same system but with smaller thermal velocity. 
Here we choose T = 0.04. For low pulling forces, we 
expect the mobility to be larger than in the case T = 1 
since the decreased thermal motion does not disturb 
the intruder travelling through the almost resting sur- 
rounding disks. For high forces, we expect the intruder 
velocity to not depend on the temperature, since in this 
case, the intruder's velocity is at least one order of mag- 
nitude higher than the thermal velocity (see Fig.[S]). 

These expectations are indeed born out by the data. 
We plot the force velocity relation for packing fractions 
77 = 0.6 and 0.775 for both temperatures in Fig. [51 
For high pulling forces F > 10'^, the intruder velocities 
for low and high temperature collapse for both packing 
fractions as expected. For rj — 0.6 and T — 0.04, the 
crossover from the linear regime to the nonlinear regime 
is shifted to smaller forces. For 77 = 0.775 and T — 0.04, 
the linear regime is not even detectable. The crossover 
is shifted by roughly two decades. To explore the linear 
regime would require forces as small as F = 10~^, at 
which the average intruder velocity would be too small 
to be detectable. 
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Fig. 6 Velocity-force relation for strong pulling forces at dif- 
ferent temperatures and e = 0.9. 



5 Conclusion and Outlook 

We have generalised an event-driven algorithm to in- 
clude friction. As a first application, we have studied the 
dynamics of an intruder pulled by an external force. In 
contrast to previous simulations [S], we consider a flu- 
idized granular medium, which is expected to undergo 
a glass transition at a packing fraction below random 
close packing [5] . We do indeed find a dramatic decrease 
of the mobility around 77 = 0.8, in agreement with pre- 
vious simulations without pulling force and consistent 
with a glass transition. For large pulling force, the data 
can be collapsed by scaling, following a power law de- 
pendence vi oc with (3 = 0.55. 

In the frictionless case the pulling force generates a 
momentum wave propagating through the sample and 
thereby compactifying it. We plan to investigate mo- 
mentum transport in more detail in the future. Fur- 
thermore the generalised event driven algorithm will be 
useful more generally in the context of frictional gran- 
ular matter fluidized by air or water flow. Work along 
these lines is in progress. 
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